#####################################################################
##Replication Code for ##############################################
##The Politics of Promotion in China's Foreign Policy Bureaucracy ###
##Tyler Jost and Yucong Li ##########################################
##The China Quarterly ###############################################
##(3) Promotion Analysis ############################################
#####################################################################


#########################################
####Load Data############################
#########################################

rm(list = ls())
load("ts_minister.RData")
load("ts_minister_connections.RData")
load("ts_vminister.RData")
load("ts_vminister_connections.RData")
load("mfa_personnel.RData")
load("mfa_appointments.RData")

crse <- function(x, df){
  d.cse <- na.omit(df[, c("id", all.vars(formula(x)))])
  fc <- d.cse$id
  m <- length(unique(fc))
  k <- length(coef(x))
  u <- estfun(x)
  u.clust <- matrix(NA, nrow=m, ncol=k) 
  for(j in 1:k){
    u.clust[,j] <- tapply(u[,j], fc, sum)
  }
  bread <-vcov(x)
  c.vcov <- bread %*% ((m / (m-1)) * t(u.clust) %*% (u.clust)) %*% bread
  result <- coeftest(x, c.vcov)
  return(result)
}

crse2 <- function(x, df){
  d.cse <- na.omit(df[, c("year", all.vars(formula(x)))])
  fc <- d.cse$year
  m <- length(unique(fc))
  k <- length(coef(x))
  u <- estfun(x)
  u.clust <- matrix(NA, nrow=m, ncol=k) 
  for(j in 1:k){
    u.clust[,j] <- tapply(u[,j], fc, sum)
  }
  bread <-vcov(x)
  c.vcov <- bread %*% ((m / (m-1)) * t(u.clust) %*% (u.clust)) %*% bread
  result <- coeftest(x, c.vcov)
  return(result)
}

##############################################
##Promotion Analysis##########################
##############################################
##Vice Minister Level#########################
##############################################

##Table 3
m1 <- glm(vm_promote ~ seniorshareabroad + as.factor(pc) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster1 <- crse(m1, vmin)
cluster1
pseudoR2.1 <- round(pR2(m1)[[4]],2)

m2 <- glm(vm_promote ~ totalmidcount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster2 <- crse(m2, vmin)
cluster2
pseudoR2.2 <- round(pR2(m2)[[4]],2)

m3 <- glm(vm_promote ~ seniorshareabroad + totalmidcount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster3 <- crse(m3, vmin)
cluster3
pseudoR2.3 <- round(pR2(m3)[[4]],2)

m4 <- glm(vm_promote ~ seniorshareabroad + totalmidcount + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
          data=vmin, family="binomial")
cluster4 <- crse(m4, vmin)
cluster4
pseudoR2.4 <- round(pR2(m4)[[4]],2)

m5 <- glm(vm_promote ~ seniorshareabroad + totalmidcount + totaltreaty + as.factor(pc) + as.factor(count_bin), 
          data=vmin[vmin$year>=1982,], family="binomial")
cluster5 <- crse(m5, vmin[vmin$year>=1982,])
cluster5
pseudoR2.5 <- round(pR2(m5)[[4]],2)

m6 <- glm(vm_promote ~ seniorshareabroad + totalmidcount + totaltreaty + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
          data=vmin[vmin$year>=1982,], family="binomial")
cluster6 <- crse(m6, vmin[vmin$year>=1982,])
cluster6
pseudoR2.6 <- round(pR2(m6)[[4]],2)

stargazer(m1, m2, m3, m4, m5, m6, type = "latex", 
          se = list(cluster1[,"Std. Error"], cluster2[,"Std. Error"], cluster3[,"Std. Error"], cluster4[,"Std. Error"],
                    cluster5[,"Std. Error"], cluster6[,"Std. Error"]),
          omit=c("Constant", "pc", "count_bin"),
          omit.stat=c("f", "ser",  "aic", "ll"),
          label=c("tab: vm_promotion"),
          covariate.labels = c("Share of Time Abroad", "International Disputes (count)", "Diplomatic Treaties (count)",
                               "Junior Priority Experience (number of posts)", "Male", "Military Background", "Started in ILD", "Higher Civilian Education", "Princeling"),
          add.lines=list(c("Party Congress Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Experience Count Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Clusters", rep(length(unique(vmin$cname)),3), 
                           length(unique(vmin$cname[is.na(vmin$juniorpriorityposts)==F])),
                           length(unique(vmin$cname[vmin$year>=1982])),
                           length(unique(vmin$cname[is.na(vmin$juniorpriorityposts)==F & vmin$year>=1982]))),
                         c("McFadden Pseudo-R$^2$", pseudoR2.1, pseudoR2.2, pseudoR2.3, pseudoR2.4, pseudoR2.5, pseudoR2.6)),
          dep.var.labels = c("Promotion To Vice Minister Rank"),
          title = "Postings, Performance, and Promotion To Vice Minister Rank",
          notes.append=FALSE, no.space=T,
          notes = c("Robust standard errors are clustered by individual.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))

##Calculate Percent Change in Odds
round(exp(coef(m1)["seniorshareabroad"]) - 1,2)*100 #87% decrease
round((exp(coef(m1)["seniorshareabroad"] * sd(vmin$seniorshareabroad, na.rm = TRUE)) - 1) * 100, 2) #48% decrease with 1 sd increase
round(exp(coef(m2)["totaltreaty"]) - 1,2)*100 #6% decrease

##############################################
##Promotion Analysis##########################
##############################################
##Minister Level##############################
##############################################

##Calculate number of individuals in the data
length(unique(minister$cname)) #139 individuals

##Table 4 Calculations
zhou.abroad <- d[d$cname=="周恩来" & d$location!="Beijing" & d$start<1949 & is.na(d$start)==F,]
zhou.home <- d[d$cname=="周恩来" & d$location=="Beijing" & d$start<1949 & is.na(d$start)==F,]

chen.abroad <- d[d$cname=="陈毅" & d$location!="Beijing" & d$start<1958 & is.na(d$start)==F,]
chen.home <- d[d$cname=="陈毅" & d$location=="Beijing" & d$start<1958 & is.na(d$start)==F,]

ji.abroad <- d[d$cname=="姬鹏飞" & d$location!="Beijing" & d$start<1972 & is.na(d$start)==F,]
ji.home <- d[d$cname=="姬鹏飞" & d$location=="Beijing" & d$start<1972 & is.na(d$start)==F,]

qiao.abroad <- d[d$cname=="乔冠华" & d$location!="Beijing" & d$start<1974 & is.na(d$start)==F,]
qiao.home <- d[d$cname=="乔冠华" & d$location=="Beijing" & d$start<1974 & is.na(d$start)==F,]

huang.abroad <- d[d$cname=="黄华" & d$location!="Beijing" & d$start<1976 & is.na(d$start)==F,]
huang.home <- d[d$cname=="黄华" & d$location=="Beijing" & d$start<1976 & is.na(d$start)==F,]

wu.abroad <- d[d$cname=="吴学谦" & d$location!="Beijing" & d$start<1982 & is.na(d$start)==F,]
wu.home <- d[d$cname=="吴学谦" & d$location=="Beijing" & d$start<1982 & is.na(d$start)==F,]

qian.abroad <- d[d$cname=="钱其琛" & d$location!="Beijing" & d$start<1988 & is.na(d$start)==F,]
qian.home <- d[d$cname=="钱其琛" & d$location=="Beijing" & d$start<1988 & is.na(d$start)==F,]

tang.abroad <- d[d$cname=="唐家璇" & d$location!="Beijing" & d$start<1998 & is.na(d$start)==F,]
tang.home <- d[d$cname=="唐家璇" & d$location=="Beijing" & d$start<1998 & is.na(d$start)==F,]

li.abroad <- d[d$cname=="李肇星" & d$location!="Beijing" & d$start<2003 & is.na(d$start)==F,]
li.home <- d[d$cname=="李肇星" & d$location=="Beijing" & d$start<2003 & is.na(d$start)==F,]

yang.abroad <- d[d$cname=="杨洁篪" & d$location!="Beijing" & d$start<2007 & is.na(d$start)==F,]
yang.home <- d[d$cname=="杨洁篪" & d$location=="Beijing" & d$start<2007 & is.na(d$start)==F,]

wang.abroad <- d[d$cname=="王毅" & d$location!="Beijing" & d$start<2013 & is.na(d$start)==F,]
wang.home <- d[d$cname=="王毅" & d$location=="Beijing" & d$start<2013 & is.na(d$start)==F,]

qin.abroad <- d[d$cname=="秦刚" & d$location!="Beijing" & d$start<2022 & is.na(d$start)==F,]
qin.home <- d[d$cname=="秦刚" & d$location=="Beijing" & d$start<2022 & is.na(d$start)==F,]
 
zhou.abroad.total <- 0
zhou.home.total <- 0

chen.abroad.total <- 0
chen.home.total <- 0

ji.abroad.total <- (1955-1950+1)
ji.home.total <- (1972-1955+1)

qiao.abroad.total <- 0
qiao.home.total <- (1974-1949+1)

huang.abroad.total <- (1965-1960+1) + (1969-1966+1) + (1971-1971+1) + (1976-1971+1)
huang.home.total <- (1960-1954+1)

wu.abroad.total <- 0
wu.home.total <- (1982-1958+1)

qian.abroad.total <- (1963-1955+1) + (1974-1972+1) + (1976-1974+1)
qian.home.total <- (1988-1976+1)

tang.abroad.total <- (1983-1978+1) + (1991-1988+1)
tang.home.total <- (1969-1964+1) + (1988-1985+1) + (2000-1991+1)

li.abroad.total <- (1977-1970+1) + (1985-1983+1) + (1995-1993+1) + (2001-1998+1)
li.home.total <- (1983-1977+1) + (1993-1985+1) + (1998-1995+1) + (2003-2001+1)

yang.abroad.total <- (1987-1983+1) + (1995-1993+1) + (2004-2001+1)
yang.home.total <- (1983-1975+1) + (1993-1987+1) + (2001-1995+1) + (2007-2004+1)

wang.abroad.total <- (1994-1989+1) + (2007-2004+1)
wang.home.total <- (1989-1982+1) + (2004-1994+1) + (2013-2007+1)

qin.abroad.total <- (1999-1995+1) + (2005-2002+1) + (2011-2010+1) + (2022-2021+1)
qin.home.total <- (1995-1992+1) + (2002-1999+1) + (2010-2005+1) + (2021-2011+1)

##Calculate averages
round(mean(c(zhou.abroad.total, chen.abroad.total, ji.abroad.total, qiao.abroad.total, huang.abroad.total, wu.abroad.total, qian.abroad.total,
             tang.abroad.total, li.abroad.total, yang.abroad.total, wang.abroad.total, qin.abroad.total)),0) #8 years
round(mean(c(zhou.home.total, chen.home.total, ji.home.total, qiao.home.total, huang.home.total, wu.home.total, qian.home.total,
             tang.home.total, li.home.total, yang.home.total, wang.home.total, qin.home.total)),0) #18 years

##Table 5
m1 <- glm(high_promote ~ seniorshareabroad + as.factor(count_bin), 
          data=minister, family="binomial")
cluster1 <- crse2(m1, minister)
cluster1
pseudoR2.1 <- round(pR2(m1)[[4]],2)

m2 <- glm(high_promote ~ allmid + alltreaties + as.factor(count_bin), 
          data=minister, family="binomial")
cluster2 <- crse2(m2, minister)
cluster2
pseudoR2.2 <- round(pR2(m2)[[4]],2)

m3 <- glm(high_promote ~ seniorshareabroad + allmid + alltreaties + as.factor(count_bin), 
          data=minister, family="binomial")
cluster3 <- crse2(m3, minister)
cluster3
pseudoR2.3 <- round(pR2(m3)[[4]],2)

m4 <- glm(high_promote ~ seniorshareabroad + allmid + alltreaties + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(count_bin), 
          data=minister, family="binomial")
cluster4 <- crse2(m4, minister)
cluster4
pseudoR2.4 <- round(pR2(m4)[[4]],2)

m5 <- glm(high_promote ~ seniorshareabroad + allmid + alltreaties + as.factor(count_bin), 
          data=minister[minister$year>=1982,], family="binomial")
cluster5 <- crse2(m5, minister[minister$year>=1982,])
cluster5
pseudoR2.5 <- round(pR2(m5)[[4]],2)

m6 <- glm(high_promote ~ seniorshareabroad + allmid + alltreaties + juniorpriorityposts + male + military + ildstart + civilcollege + princeling + as.factor(count_bin), 
          data=minister[minister$year>=1982,], family="binomial")
cluster6 <- crse2(m6, minister[minister$year>=1982,])
cluster6
pseudoR2.6 <- round(pR2(m6)[[4]],2)

stargazer(m1, m2, m3, m4, m5, m6, type = "latex", 
          se = list(cluster1[,"Std. Error"], cluster2[,"Std. Error"], cluster3[,"Std. Error"],
                    cluster4[,"Std. Error"], cluster5[,"Std. Error"], cluster6[,"Std. Error"]),
          omit=c("Constant", "count_bin"),
          omit.stat=c("f", "ser", "aic", "ll"),
          label=c("tab: min_promotion"),
          covariate.labels = c("Share of Time Abroad", "International Disputes (count)", "Diplomatic Treaties (count)", "Junior Priority Experience (number of posts)", "Male", "Military Background", "Started in ILD", "Higher Civilian Education", "Princeling"),
          add.lines=list(c("Experience Count Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Clusters", rep(length(unique(minister$year)),3),
                           length(unique(minister$year[is.na(minister$juniorpriorityposts)==F])),
                           length(unique(minister$year[minister$year>=1982])),
                           length(unique(minister$year[is.na(minister$juniorpriorityposts)==F & minister$year>=1982]))),
                         c("McFadden Pseudo-R$^2$", pseudoR2.1, pseudoR2.2, pseudoR2.3, pseudoR2.4, pseudoR2.5, pseudoR2.6)),
          dep.var.labels = c("Appointment to Minister Rank"),
          title = "Postings, Performance, and Appointment to Minister Rank",
          notes.append=FALSE, no.space=T,
          notes = c("Robust standard errors are clustered by year.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))

##############################################
##Network Analysis############################
##############################################

# Descriptive Figures
round(mean(vmin_connections$totalconnections_count, na.rm = T),0) #19 connections
vmin_connections$totalconnections_count[vmin_connections$cname=="唐家璇" & vmin_connections$year==1993]  #45 connections
round(mean(vmin_connections$totalconnections_count[vmin_connections$year==1993]),0) #14 connections

# Vice Minister Level
connections1 <- glm(totalconnections_count ~ seniorshareabroad + male + military + civilcollege + princeling + as.factor(pc) + as.factor(count_bin), 
                    data=vmin_connections)
crse(connections1, vmin_connections)
round(coef(connections1)["seniorshareabroad"]) #14% decrease

# Minister Level
connections2 <- glm(totalconnections_count ~ seniorshareabroad + male + military + civilcollege + princeling + as.factor(count_bin) + as.factor(year), 
                    data=minister_connections)
crse(connections2, minister_connections)
round(coef(connections2)["seniorshareabroad"]) #7% decrease

